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Abstract 

We compute masses of D mesons, Dg mesons and charmonium states using Wilson twisted 
mass lattice QGD. We present results for spin J = 0,1,2,3, parity V = —, + and in case of 
charmonium also charge conjugation C = —Computations are based on quark-antiquark 
creation operators and performed at three different unphysically heavy u/d quark masses allow¬ 
ing an extrapolation to the physical u/d quark mass. Within combined statistical and systematic 
errors, which are around 2%... 3%, our results agree with available experimental results. Par¬ 
ticular focus is put on the = 1+ mesons Hi (2430) and Hi (2420). We separate and classify 
these close-by states according to the total angular momentum of their light degrees of freedom, 
j « 1/2 and j « 3/2. This is a first important step to study decays B^*'^ Di + l + v, for which 
a long-standing conflict between theory and experiment exists. 




1 Introduction 


More than fifty D meson, Dg meson and charmonium states have been observed in experiments 
[1]. Several of them seem to be theoretically well-understood, e.g. the = 0“ pseudoscalar 
and = 1“ vector ground states. There are, however, open questions regarding some of the 
more recently found excitations. A prominent example is the charmonium-like state Al(3872) 
first observed by Belle [2] , which is most likely not an ordinary quark-antiquark state and whose 
structure is under debate. The situation is similar for several other charmonium-like X states. In 
particular the electrically charged states, e.g. X(3900)^ or X(4020)^, cannot be simple cc pairs. 
One rather assumes a four quark structure, e.g. a mesonic molecule or a diquark-antidiquark pair. 
The situation for open charm mesons is related, even though fewer D and Dg states have been 
observed. For example the positive parity mesons T>*q(2317) and L>si(2460) first reported by 
BaBar [3] and CLEO [3], respectively, are unexpectedly light. Again this could be an indication 
that these states are not just quark-antiquark pairs, but are composed of two quarks and two 
antiquarks, a scenario at the moment neither established nor ruled out. 

There are many interesting and to some extent successful approaches to study D mesons, Dg 
mesons and charmonium states theoretically, e.g. quark models [5], effective theories respect¬ 
ing QCD symmetries [6] or Dyson-Schwinger and Bethe-Salpeter equations [7], to just name a 
few. Of course, it would be highly desirable to understand these mesons and their properties 
starting from first principles, i.e. the QCD Lagrangian, without any assumptions, model sim¬ 
plifications or truncations. The corresponding and commonly used method to achieve that goal 
is lattice QCD, a numerical method to compute QCD observables. In principle, it allows to 
investigate and to quantify all possible sources of systematic error. Lattice meson spectroscopy 
is, however, a challenging task, where many problems have currently only partly been solved. 
On the one hand there are issues concerning lattice QCD in general. For example simulations 
with physically light u/d quarks are extremely demanding with respect to high performance 
computing resources. Similarly, to remove discretization errors one has to study the continuum 
limit, which requires simulations at several different lattice spacings, again a very challenging 
task with respect to computational resources. On the other hand there are problems specific 
to lattice hadron spectroscopy. An example is the investigation of states, which can decay into 
lighter multi-particle states. Such states should theoretically be treated as resonances and not 
as stable quark-antiquark states, which is technically extremely difficult, even for simple cases, 
where only a single decay channel exists. Examples are Dq( 2400) and Di(2430) with quan¬ 
tum numbers = 0"*“ and = 1"*“. Similarly, it is very challenging to study mesons, which 
might have a structure more complicated than a simple quark-antiquark pair, e.g. candidates 
for tetraquarks or hybrid mesons. While there has been a lot of impressive progress regarding 
lattice hadron spectroscopy within the last couple of years, there is certainly still a lot of room 
for improvement. Simple states, in particular pseudoscalar and vector ground state mesons, 
have, meanwhile, been studied very accurately, including simulations at or extrapolations to 
physically light u/d quark masses and the continuum limit. On the other hand, the majority 
of studies concerned with parity, radial and orbital excitations are still at a more exploratory 
stage, i.e. have quite often been performed at unphysically heavy quark masses or at a single 
finite lattice spacing. Recent reviews discussing the status of lattice QCD computations of D 
and Dg mesons and of charmonium are [8] and mm, respectively. 

The most common approach to compute meson masses using lattice QCD is to employ meson 
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creation operators, which are composed of a quark and an antiquark operator, and to extract 
meson masses from the exponential decay of corresponding correlation functionqj. This strategy 
yields accurate and solid results for mesons, which resemble quark-antiquark pairs and which 
are quite stable, i.e. many of the low-lying states in the D meson, Dg meson and charmonium 
sector. Recent lattice QCD papers following this strategy to compute masses and spectra of D 
and Dg mesons and of charmonium are [IH |T3l [HI |T5l [T6l [JTl jlHl [191 [IS [HI Ell [23l EH EH 
I26j . Rigorous treatments of more complicated mesonic systems like the previously mentioned 
unstable Dq( 2400) and Hi (2430) mesons or any of the tetraquark candidates, e.g. X(3900)^ 
or X(4020)^, require more advanced techniques, in particular the implementation of meson 
creation operators composed of two quark and two antiquark operators and possibly studies of 
the volume dependence of the masses of corresponding scattering stated. Examples of recent 
lattice papers exploring and using such techniques to study specific D, Dg or charmonium states 
are [281 EH EQl El EH EH El EH EHEZ] . 

The goal of this paper is to compute the masses of several low-lying D meson, Dg meson and 
charmonium states using Wilson twisted mass lattice QCD with 2-|-l-|-l dynamical quark fla¬ 
vors. One of the main advantages of this particular discretization of QCD is automatic 0{a) 
improvement, i.e. discretization errors appear only quadratically in the small lattice spacing a 
and are, hence, quite small. From a technical point of view we employ a large variety of quark- 
antiquark meson creation operators and are, hence, able to study total angular momentum 
J = 0,1, 2,3, parity V = —, + and in case of charmonium charge conjugation C = Com¬ 

putations are performed at several unphysically heavy u/d quark masses (corresponding pion 
masses m,r ~ 276 MeV, 315 MeV, 443 MeV), which allow extrapolations to the physical point. 
At the moment computations are, however, restricted to a single lattice spacing a ~ 0.0885 fm, 
i.e. we are currently not able to perform a continuum extrapolation. Nevertheless, by using 
two different Wilson twisted mass discretizations of the meson creation operators we are able 
to crudely estimate the magnitude of discretization errors associated with our resulting meson 
masses. Computations at smaller lattice spacings and corresponding continuum extrapolations 
are planned for the near future and will be part of an upcoming publication. 

As mentioned above some D meson, Dg meson and charmonium states are quite unstable or 
might have a structure much different from a quark-antiquark pair. Even though we present 
results for these states in the following, a rigorous treatment might require more advanced 
techniques, in particular the inclusion of four-quark creation operators as discussed above. We 
are in the process of developing such techniques using a similar lattice QCD setup [38l |39l SO] • 
The techniques and results presented in this paper are an important prerequisite for such more 
advanced computations. Our long-term goal is the computation of the low-lying D and Dg meson 
and charmonium spectra as fully as possible with all sources of systematic error removed or 
quantified (in particular computations at physically light u/d quark masses including continuum 
extrapolations), using quark-antiquark creation operators supplemented, whenever necessary, by 
four-quark creation operators. 

In this work we also study the structure of the two lightest D mesons with = I"*", Hi(2430) 
and Hi (2420). These states are similar in mass, but their structure is quite different. One of 
them has j ~ 1/2, while the other has j ~ 3/2, where j denotes the total angular momentum 
of the light quark and of gluons. We demonstrate, how to resolve and classify both states 

^For a basic introduction to lattice hadron spectroscopy cf. m- 

^For a basic introduction on how to study resonances using lattice QCD cf. m 
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from computations based on a single = 1^ correlation matrix. The eigenvector components 
obtained by solving a generalized eigenvalue problem provide linear combinations of meson 
creation operators suited to specifically excite T)i(2430) and Di(24:20). This is an important 
hrst step to study decays B^*') —>• D**+l+u, for which there is a persistent conflict between theory 
and experiment regarding branching ratios (cf. |41] for a detailed discussion). Our techniques 
and results can be used to extend recent lattice QCD computations of these decays, where D** 
has been restricted to = 0''',2+, but did not include T)i(2430) and i4i(2420) [4^ 143] . 

Parts of this work have been presented at recent conferences [441I4S1I45]. 

This paper is structured as follows. In section [2] we introduce the 2+1+1 flavor Wilson twisted 
mass lattice setup. In section [3| and section |4| we discuss the technical aspects of our com¬ 
putations, in particular the employed meson creation operators, their corresponding quantum 
numbers and how we compute correlation matrices and meson masses. In section [5] we present 
our results, first for D and Dg mesons in section EH then for charmonium in section [+21 These 
results are summarized in plots and tables in section El where we also give a brief outlook. 


2 Lattice QCD setup 

2.1 Gauge link ensembles, sea quarks 

We use gauge link configurations generated with 2+1+1 dynamical quark flavors by the Euro¬ 
pean Twisted Mass Collaboration (ETMC) [471 SSI SHI ED]. The gluonic action is the Iwasaki 
gauge action |5l]. For the light degenerate {u,d) quark doublet the standard Wilson twisted 
mass action 

*S'iight[x^'\x^'\C = +in'y5T3^x^^\x) (1) 


has been used |52|, for the heavy (c, s) sea quark doublet the Wilson twisted mass formulation 
for non-degenerate quarks 

5'heavy[x^^\x^^\t^] = (^Dw{mo) + (2) 


m- Hw denotes the Wilson Dirac operator, 

Dw(mo) = ^(7 m(v^ + V*) - +mo, (3) 

x(0 = ^{h) ^ 

are the quark fields in the so-called twisted basis and ri 
and T 3 denote the first and third Pauli matrix acting in flavor space. At maximal twist physical 
quantities, e.g. meson masses, are automatically 0(a} improved [54j- The tuning has been done 
by adjusting mo such that the PCAC quark mass in the light quark sector vanishes (cf. [4H| for 
details). For a review on Wilson twisted mass lattice QCD we refer to [55] . 

In this work we use three ensembles with gauge coupling /3 = 1.90, which amounts to a lattice 
spacing o « 0.0885fm (scale setting via the pion mass and the pion decay constant [HD])- The 
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ensembles differ in the unphysically heavy u/d quark mass /x = 0.0030, 0.0040, 0.0080 corre¬ 
sponding to tUtt ~ 276 MeV , 315 MeV , 443 MeV. The s and the c quark masses are represented 
by /Xct = 0.150 and /x^ = 0.190. These values have been chosen such that the lattice QCD results 
for and for quantities, which depend only weakly on the light u/d quark mass, 

are close to the corresponding physical values [121153 ES] • Details of these gauge link ensembles 
are collected in Table [3 


ensemble 


{L/af X T/a 




a 

(fm) 

W-TT 

(MeV) 

#of 

configurations 

A30.32 

1.90 

32^ X 64 

0.0030 

0.150 

0.190 

0.0885 

276 

1200 

A40.32 


32^ X 64 

0.0040 




315 

800 

A80.24 


24^ X 48 

0.0080 




443 

1700 


Table 1: Gauge link ensembles {{L/af’ x T/a: number of lattice sites; # of configurations: 
number of gauge link configurations used). 


2.2 Valence quarks 

For the light degenerate (xx, d) valence quark doublet we use the same action, which was used to 
simulate the corresponding sea quarks, i.e. the action ([T|). 

For the heavy s and c valence quarks we use twisted mass doublets of degenerate quarks, i.e. 
a different discretization as for the corresponding sea quarks. We use the action ([T]) with the 
replacements = (x^^^\x^* ^)! and —>■ 

respectively. We do this, to avoid mixing of strange and charm quarks, which inevitably 
takes place in a unitary non-degenerate Wilson twisted mass setup, and which is particularly 
problematic for observables containing charm quarks, e.g. masses of D and Dg mesons and of 
charmonium (cf. [53 I5B] for a detailed discussion of these problems). 

The degenerate valence doublets allow two realizations for strange as well as for charm quarks, 
either with a twisted mass term +ifis,c 75 (he. ^ or ^) or —ifJ,s,c 75 (i-e. x^* ^ or x^^ ^)- 
For a quark-antiquark meson creation operator, e.g. X^^HsX^^^ the sign combinations (-1-, —) 
and (—,-|-) for the antiquark and the quark are related by symmetry, i.e. the corre¬ 
sponding correlation functions are identical. These correlation functions differ, however, from 
their counterparts with sign combinations (-|-,-|-) and (—,—) due to different discretization er¬ 
rors. In section|5]we will show for each computed meson mass both the (-I-, —) = (—, -|-) and the 
(-I-, -k) = (—, —) result. The mass differences are O(a^), because of automatic 0{a) improvement 
at maximal twist. These differences provide a first estimate of the magnitude of discretization 
errors at our currently used lattice spacing a ^ 0.0885 fm. 

The tuning of the valence quark masses /x^ and fic is discussed in section 14.31 


4 












3 Meson creation operators and trial states 

3.1 Meson creation operators in the continnum 

In the continuum a quark-antiquark operator creating a mesonic trial state with definite quantum 
numbers (total angular momentum J, parity V, charge conjugation C), when applied to the 
vacuum | 0 ), is 

{l/\/V) f d^r projects to vanishing total momentum {V is the spatial volume), i.e. realizes a 
meson at rest. To allow for orbital angular momentum between the antiquark and the quark 
they are spatially separated. i|Ar|=i?denotes an integration over a sphere of radius 
R, which is the distance between the antiquark and the quark. r(Ar) is a suitable combination 
of spherical harmonics and 7 matrices (cf. Tabled column ‘T(n), pb”), which combines orbital 
angular momentum and the two quark spins to total angular momentum J and determines parity 
V and in case of identical quark flavors charge conjugation C. [/(r; r + Ar) is a straight gluonic 
parallel transporter connecting the antiquark and the quark in a gauge invariant way. For D 
mesons e.g. = uc, for Dg mesons e.g. = sc and for charmonium = cc. 


d^Ar U{r;r + Ar)r(Ar)?/)^^^(r -|- Ar). (4) 


3.2 Meson creation operators in Wilson twisted mass lattice QCD 

Our lattice meson creation operators are of similar form, 

^ ^ t/(n;n + An)r(An)xl=>(n + An), (5) 

V n An=±ea;,±ej,,±e2 

where the integration over a sphere with center at r has been replaced by the sum over the 
six neighboring lattice sites of n and U (n; n + An) denotes the link between n and n + An. 
Moreover, physical basis quark operators have been replaced by their twisted basis 

counterparts 


3.2.1 Physical basis and twisted basis 

In the continuum the relation between the physical and the twisted basis is given by the twist 
rotation 

— exp ^i75r3a;/2^ j exp ^i75r3a;/2^ (6) 

with the twist angle ui, where oj = '7r/2 at maximal twist, x^'^^ denotes either the light 
doublet x^^^ = the strange doublet x^^^ = (x^*^\x^* ^) or the charm doublet 

X*''^^ = X^^ ^) (cf- also section [ 22 ]) . 

When transforming a twisted basis quark bilinear x^^^hx^^^ as e.g. appearing in (|5|) to the 
physical basis or vice versa, the result depends not only on F, but also on the flavor combination, 
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continuum 

twisted mass lattice QCD 

index 

r(n), pb 

J 

VC 

tb, (+,+) 

tb, (±,±) 

0^ ® 0^ - 


1 

75 



pb 

±*75 X 



2 

3 

7075 

1 


■ 

±i75X 

pb 

pb 

±*75 X 

Ai (8)yli 


4 

70 

0 

■ 

+Z75X 

pb 


Ai 

5 

757jnj 

— 

±*75 X 

pb 


6 

7o757inj 


-+ 

pb 

±*75 X 

Ti (gjTi 


7 

71%- 


++ 

±*75 X 

pb 


8 

707jnj 


++ 

pb 

±*75 X 



1 

71 


B 

±*75 X 

pb 



2 

3 

7o7i 

7571 


■ 

pb 

±*75 X 

±*75 X 

pb 

■ 


4 

707571 


■ 

pb 

±*75 X 



5 

ni 


B 

pb 

±*75 X 



6 

7 

7oni 

75ni 


■ 

±*75 X 
pb 

pb 

±*75 X 

Ai(E)Ti 


8 

7o75ni 

1 

19 

±*75 X 

pb 


Ti 

9 

(n X 7)1 

++ 

±*75 X 

pb 


10 

7o(n X 7)1 


++ 

pb 

±*75 X 

Ti(E)Ti 


11 

75(n X 7)1 


— 

±*75 X 

pb 


12 

7075 (n X 7)1 


-+ 

pb 

±*75 X 



13 

7i(2n^ -112 -n^) 


— 

±*75 X 

pb 



14 

7071 (2nf - n| - n|) 


— 

pb 

±*75 X 

Ti0E 


15 

757i(2nf -n| -n2) 


++ 

±*75 X 

pb 


16 

707571 (2nf - n| - n|) 


+- 

pb 

±*75 X 



1 

(n^ + ni- 2ni) 


++ 

pb 

±*75 X 



2 

7o(nf + n2 - 2n|) 


+- 

±*75 X 

pb 

Ai0E 


3 

75(nf +n| - 2n|) 


-+ 

pb 

±*75 X 


4 

7075(nf + nl- 2n|) 

9 

-+ 

±*75 X 

pb 


E 

5 

(7ini + 72112 - 273113) 


91 

±*75 X 

pb 


6 

70 (71111 + 72112 - 273113) 


19 

pb 

±*75 X 

Ti(E)Ti 


7 

75(71111 + 72112 - 273113) 


H 

±*75 X 

pb 


8 

7075(71111 + 72112 - 273113) 



pb 

±*75 X 



1 

(73112 + 72113) 


++ 

±*75 X 

pb 



2 

70(73112 + 72113) 


++ 

pb 

±*75 X 

Ti (gjTi 


3 

75(73112 + 72113) 


— 

±*75 X 

pb 


4 

7075(73112+72113) 

9 

-+ 

pb 

±*75 X 


T 2 

5 

71(112-113) 


— 

±*75 X 

pb 


6 

7 

7071 (ni - ni) 

7571 (ni - ni) 


++ 

pb 

±*75 X 

±*75 X 

pb 

Ti^E 


8 

707571 (ni - ni) 


+- 

pb 

±*75 X 




Table 2: Meson creation operators. 
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i.e. whether and are upper components (twisted mass term +i/i 75 ) or lower components 
(twisted mass term of twisted basis doublets. For example 

^ = x^“^exp(^ + i75a;/2^75exp(^-i75u;/2^X^'' ^ \ (7) 

while 

exp + ij 5 U}/ 2 ^ 75 exp + ij 5 U}/ 2 j = 

= X^“^exp (^ + Z75a;^75x^''’^^ X^'^H+ijbhbX^"^'' = (8) 

In the columns “tb, (=h, ^)” and “tb, (zh, zb)” of Table[2]we list for all flavor combinations (+ and 
— denote the signs in front of the twisted mass terms for x*-^^ and x*'^^) and all F combinations 
of our meson creation operators, how physical and twisted basis are related, “pb” indicates that 
the twisted basis F is the same as the physical basis F (cf. e.g. ([7])), while “zbzysx” denotes that 
the physical basis F has to be multiplied from the left with zbiys to obtain the corresponding 
twisted F (cf. e.g. ([ 8 |)). 

At finite lattice spacing the twist rotation only holds for renormalized operators, i.e. for bare 
lattice quark operators the twist rotation ([ 6 ]) is only an approximate relation. Nevertheless, it is 
possible to unambiguously interpret states obtained from correlation functions of twisted basis 
meson creation operators in terms of QCD quantum numbers as we will explain and demonstrate 
in section Ol 


3.2.2 Isospin, parity and charge conjngation 

Isospin I and parity V are symmetries of QCD. While in Wilson twisted mass lattice QCD the 
2 ; component of isospin is still a quantum number, I and V are broken by 0(a) due to the 
Wilson term appearing in the twisted mass actions ([I]) and ([3]). Only a 

specific combination of both symmetries, light flavor exchange u d combined with parity, is 
still a symmetry. We denote this symmetry by acting on the light twisted basis quark 

doublet x*'^^ = according to = 7 oTix^^^- Similarly, = yonx^^^ 

and X^) = yoTix^'^^- Note that 7 ^ 0. In general, it is, therefore, not possible to 

classify states according to and at the same time. 

For D mesons we use trial states ( 2 )|with defined e.g. is 

suited for D mesons with = +1/2. There are eight appropriate flavor combinations for D 
mesons, where the four with opposite signs in front of the twisted mass terms (denoted by (±, +) 
throughout the paper), 

yWy(c+) yWyi''") yi''~iyi“i (9) 

A A ^ A A ?A A ?A A? J 


are related by symmetry and yield identical correlation functions. Similarly the four flavor 
combinations with identical signs in front of the twisted mass terms (denoted by (+,+)), 






^(c+)xW 




( 10 ) 
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also yield identical correlation functions. However, at finite lattice spacing (±,^) and (±,±) 
correlation functions slightly differ. As a consequence D mesons computed on the one hand 
with a (±,±) and on the other hand with a (±,d-) flavor combination, but which are otherwise 
identical, will differ in mass. Due to automatic 0{a) improvement of Wilson twisted mass lattice 
QCD at maximal twist, this mass splitting will be proportional to a^, i.e. is expected to be rather 
small and will vanish quadratically, when approaching the continuum limit. Even though we 
consider only a single lattice spacing a ^ 0.0885 fm throughout this work, the splitting between 
(zh,^) and (±,zh) flavor combinations will provide a crude estimate of the magnitude of lattice 
discretization errors associated with the resulting meson masses. 

Since parity is not a symmetry, there is no rigorous separation between V = + and V = — 
states in Wilson twisted mass lattice QCD. Corresponding correlation functions, e.g. between the 
meson creation operators (continuum quantum numbers = 0“) and 

(continuum quantum numbers = 0^), are 0{a), but do not vanish. 

Identical considerations apply for Dg mesons, when replacing {u,d) 

For charmonium states there is an additional symmetry, charge conjugation C, where = 

7o72(x^'^^)^- For charmonium creation operators there are two appropriate flavor combinations, 

y(c+)y(c+) 

which are again related by symmetry. One might also consider flavor combinations with opposite 
signs in front of the twisted mass terms, 

y(c+)y(c-) yF'lyl^"^) 

In this case C is not a symmetry, but C combined with pd™) denoted by Copd™). Note, however, 
that a rigorous treatment of charmonium states with flavor combinations (II2p is not possible, 
because disconnected diagrams are excluded by construction. Since disconnected diagrams are 
expected to be negligible compared to statistical errors and, therefore, ignored in this work (cf. 
also the discussion in section 15.2p , computations using (fT^ are still useful to get a hrst estimate 
of lattice discretization errors (as discussed above, the differences between (±,zh) and (±,=p) 
charmonium masses are proportional to a?). 

3.2.3 Rotational symmetry and total angular momentum 

On a cubic lattice rotational symmetry is reduced to symmetry with respect to cubic rotations. 
There are only hve different irreducible representations of the cubic group O (labeled by Ai, Ti, 
E, T 2 , A 2 ), each corresponding to an inhnite number of SO(3) irreducible representations in the 
continuum (labeled by a non-negative integer referring to e.g. spin, orbital angular momentum 
or total angular momentum S, L, J = 0,1,2,...): 


Ai 


angular momenta 0,4,.. 

Ti 


angular momenta 1,3,4, 

E 


angular momenta 2,4,.. 

T 2 

—)• 

angular momenta 2,3,4, 

A 2 


angular momenta 3,... 




The two quark spins of the meson creation operators ([5|) can be combined via 7 matrices to 
S' = 0,1, which corresponds to the Ai (singlet) and the Ti (triplet) representation. The spin is 
coupled to orbital angular momentum L, where we can access with our choice of summing over 
six neighboring lattice sites in ([5]) the cubic representations Ai, Ti and E. Therefore, the total 
angular momentum O representations of our lattice meson creation operators are 

Ai® Ai®Ti®Ti®Ti®Ti®E®E®T2®T2. 

^ V 

total angular momentum J 

(13) 

The meson creation operators listed in Table [2] are sorted and organized according to these 10 
multiplets (cf. the column “O'^ ® —)• O'^”). 

3.3 Smearing of gauge links and quark fields 

To enhance the overlap of trial states |fl) to low lying meson states, we use stan¬ 

dard smearing techniques. This allows to read off meson masses from the exponential decay 
of correlation functions at rather small temporal separations, where the signal-to-noise ratio is 
favorable. 

Smearing is done in two steps. First we replace spatial gauge links by their APE smeared 
counterparts. Then we use Gaussian smearing on the quark fields which 

resorts to the APE smeared spatial links. The parameters we have chosen are AApe = 10, 
CKAPE = 0.5, iVcauss = 30 and Kcauss = 0.5. This corresponds to a Gaussian width of the 
smeared quark fields of approximately 2.7 x a ~ 0.24fm (cf. [59] for detailed equations). 

Smearing is a symmetric operation with respect to cubic rotations and spatial reflections. There¬ 
fore, it does not change the quantum numbers J and V generated by the corresponding meson 
creation operators as listed in Table [2j 


(^Ai © Ti) © (^Ai © Ti © e) 

spin S orbital angular momentum L 


4 Computation and analysis of correlation matrices 


4.1 Computation of correlation matrices 

Eor each twisted mass sector characterized by flavor the cubic representation and, 

in case of charmonium, either C (for twisted mass signs (±, ©)) or (for twisted mass 

signs (©,©)), we compute temporal correlation matrices of meson creation operators 


(*) = {fi| (s(0“?©B))) *(‘) (•S(Oj.’"|a7B,)) (0)|i!) 


(14) 


j and k label the rows and columns of a correlation matrix or, equivalently, are indices of the 
meson creation operators entering a correlation matrix (cf. Table O column “index”). 

®The quantum number associated with a twisted basis meson creation operator from Table [21 column 

“tb, (±, =p)” is the product of the V and C quantum numbers also listed in Table (21 column “PC”. 
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indicates that APE smeared gauge links and Gaussian smeared quark fields are used for the 
meson creation operators (cf. section [331) • the computations we use a generalization of the 
one-end trick, which is explained in detail in appendix (A) 

Since parity is only an approximate symmetry in twisted mass lattice QCD, we consider corre¬ 
lation matrices of meson creation operators with both V = + and V = —. 


For D and mesons: 


- ^1, E, T2 : 8 X 8 correlation matrices. 

- Ti: 16 X 16 correlation matrix. 

For charmonium and twisted mass signs (=h,zb), i.e. s-nd x^'^ ^x^^ 

- Ai, E and C = +: 6 x 6 correlation matrices; 

Ai, E and C = —: 2x2 correlation matrices. 


— Ti and C = + 
Ti and C = — 

— T 2 and C = + 
T 2 and C = — 


6 x 6 correlation matrix; 

10 X 10 correlation matrix. 

4x4 correlation matrix; 
4x4 correlation matrix. 


For charmonium and twisted mass signs (±,q:), i.e. ^ x^'^ 

- Ai, E andCoph”) = +: 4 X 4 correlation matrices; 

Ai, E and C o ph™) = —; 4 X 4 correlation matrices. 

- Ti and C o pd”) = -K: 10 X 10 correlation matrix; 

6 x 6 correlation matrix. 


Ti and C o = - 

- T 2 and C o pdm) = + 
T 2 and C o = - 


6 x 6 correlation matrix; 
2 x 2 correlation matrix. 


In Table [T] we list for each ensemble the number of gauge link configurations used for the 
computation of the correlation matrices Cp ..p^.^( 1 )^( 2 ). The four stochastic sources needed for 
the one-end trick (cf. eq. (|27l) l are located on a timeslice, which is randomly chosen for every 
gauge link conhguration. We only use a single set of four stochastic timeslice sources, i.e. a 
single sample, for each gauge link conhguration. 

We have checked the computation of the correlation matrices Cp^.p^.^( 1 )^( 2 ) (t) by testing the 
symmetries twisted mass 75 hermiticity, twisted mass parity twisted mass time reversal, 

charge conjugation C and cubic rotations on the level of correlation matrix elements. In a second 
step we have averaged the elements related by these symmetries, to improve the signal-to-noise 
ratio. 


4.2 Determination of meson masses and assignment of parity 

We determine meson masses from the correlation matrices Cp^.p^.^( 1 )^( 2 ) specihed in the previous 
subsection (in the following denoted by Cjk = simplicity). 


10 




In a first step we replace 


Cjk{t) 


Cjk{t) 


C{t) 

\/ C'jj{t — CL^Ckkii — cl) 


(15) 


This amounts to a correlation matrix C{t) with meson creation operators Oj = Oj/Cjj {t = a), 
i.e. operators, which are normalized differently than operators ([5]), but are otherwise identical. 
Such a normalization clearly does not change the exponential decay of correlation matrix ele¬ 
ments, i.e. the meson masses we are interested in. However, it corresponds to trial states Oj\i}), 
which have a similar norm. This is important both for a correct assignment of parity and for a 
meaningful interpretation of the structure of the state associated with an extracted meson mass. 

We solve generalized eigenvalue problems 

C{t)v^^\t) = X^^\t)C{to)v^^\t) (16) 

with to = a 0 (for a detailed discussion of this generalized eigenvalue problem cf. [60] and 
references therein). For an x correlation matrix C{t) one obtains N eigenvalues and 

N corresponding eigenvectors (t), n = 0, N — 1. 

From the eigenvalues we compute N effective masses m^{t) by solving 

cosh(mf(t)(r /2 -t)) ^ Afo)(t) 

cosh(m®®(t)(r/2 — (t -I- a))) A(fo(t -|- a) 

While for mn{T/2 — t) S> 1 this equation is equivalent to the commonly known definition 


mf{t) 



a 


( Afo)(t) \ 
VAfo)(t+ a)y’ 


(18) 


(HZD yields in contrast to (1181) plateau-like effective masses m®®(t) ~ const also for large temporal 
separations t in the region t ~ r/2. The low-lying meson masses in the sector defined by 
the operators of the correlation matrix C, i.e. by Fj and by then determined by 

performing uncorrelated x^-minimizing fits of constants nin to the plateau-like regions of m^{t) 
at sufficiently large t. 

As an example we consider a 2 x 2 correlation matrix with quark flavors = x^* 

(i.e. twisted mass signs (—, -|-)), representation Ai and Fq = 75 (operator index 1 in Table[2|) 
and Fi = 1 (operator index 3 in Table [2|). The two resulting effective masses are plotted 

in Figured) (top). Clearly, there are plateaus at large temporal separations t. The fitting ranges 
for the constants rUn are indicated by straight orange lines. The resulting values for and the 
corresponding x^/dof< 1 are also included in the plot. 

The eigenvectors allow to assign parity to the extracted meson masses and to make 

qualitative statements about the structure of the corresponding states. The absolute values of 

^Theoretical arguments given in [60] suggest to choose to > tl2, since then unwanted contributions of excited 
states are strongly suppressed, in particular, when using large correlation matrices. In practice we find effective 
masses, which are essentially independent of to, but with statistical errors increasing for to > a. 
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Figure 1; Determination of meson masses and assignment of parity for Dg mesons with J = 0: 
2 x 2 correlation matrix with quark flavors = x^^ \ 0 “^ representation Ai and 

Fj G { 75 , 1 } (A30.32 ensemble). 


the j-th entry of the normalized vector 






(19) 


indicates, to which extent the meson creation operator Oj associated with the j-th line and 
j-th column of the correlation matrix C{t) excites the state corresponding to the nth extracted 
meson mass m„. 

In Figured] (bottom, left) the “operator content” of the ground state (corresponding to 

the extracted mass mo) is plotted as a function of the temporal separation t. At large t the 
meson creation operator with F = 75 is clearly dominating (~ 95%), while the meson creation 
operator with F = 1 contributes on a negligible level (~ 5%). Since F = 75 corresponds to 
negative parity, we assign the quantum number V = — to the ground state (i.e. identify that 
state as the Dg meson). From a similar plot (Figured! [bottom, right]) one can infer that the 
first excitation, which is dominated by F = 1, has V = +, i.e. should correspond to D*q. 
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4.3 Extrapolating meson masses to physical strange and charm valence quark 
masses 

For each ensemble we compute rriK at two different values of the strange valence quark mass, 
and both in the region of the physical value. By means of a linear extrapolation we 
then determine iJ,s,phys such that phys agrees with the experimental result 

= 0.477 GeV^ [1]: 


^^s,phys = = 2mK{^lsf - ml ( 20 ) 

— A y^ls,2) 

(for rrixilJ-s) and m,r we use (±,4=) twisted mass sign combinations, which are known to yield 
smaller discretization errors 16I1E2]). Since in leading order chiral perturbation theory 2ml^—ml 
is independent of the light u/d quark mass, /r^^phys should be very close to the physical strange 
quark mass. Results are collected in Table [3l For ensemble A80.24 the procedure is illustrated 
in Figure[2](top left; red points correspond to ^s,i and ^s, 2 , the black dashed lines to /is^phys)- To 
verify the validity of these linear extrapolations, we performed for ensemble A80.24 additional 
computations with a third value of the strange valence quark mass, The corresponding 

result 2mK{^is)‘^ — (tke blue point in Figure [2] [top left]) is consistent with the linear 

extrapolation. 


ensemble 

X 

Ms,a: 

2ml^ — ml in GeV^ 

l^c,x 

mo in GeV 

A30.32 

1 

0.018750 

0.454(2) 

0.22700 

1.782(2) 


2 

0.022800 

0.551(2) 

0.27720 

2 .002(2) 


phys 

0.01969(8) 

0.477* 

0.2459(4) 

1.865* 

A40.32 

1 

0.018750 

0.449(1) 

0.23887 

1.841(2) 


2 

0.023220 

0.555(1) 

0.27678 

2.007(3) 


phys 

0.01994(5) 

0.477* 

0.2443(5) 

1.865* 

A80.24 

1 

0.018749 

0.452(2) 

0.22999 

1.820(3) 


2 

0.023280 

0.561(2) 

0.29299 

2.080(2) 


3 

0.016844 

0.406(2) 

0.21438 

1.748(2) 


phys 

0.01979(8) 

0.477* 

0.2408(5) 

1.865* 


Table 3; Determining physical strange and charm valence quark masses (*: experimental results 
for 2m^o — ml^ and m do from [T] ). 


Proceeding in an analogous way we determine a charm valence quark mass /ic,phys for each 
ensemble, which is close to the physical charm quark mass, this time using niD- 


Mc,phys 


/^c,2 (/^c,l /^c,2) 


»^D,phys — mD{fic,2) 

moifj-c,!) - mDilJ-c,2) 


( 21 ) 


with mD,phys = mDO = 1.865GeV [T] (cf. Tableland Figured [top right]). Again we tested 
the quality of the linear extrapolation for ensemble A80.24 with a third value of the charm 
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determination of the physicai strange valence quark mass 


determination of the physical charm valence quark mass 



t^s t^c 


Dg mesons, twisted mass sign combination (+,-) Dg mesons, twisted mass sign combination (+,-) 





M-c 


Figure 2: (top): Determining physical strange and charm valence quark masses by means of 
linear extrapolations (A80.24 ensemble), (bottom): Linear extrapolations of various Dg meson 
masses both in and to physical strange and charm valence quark masses. 


valence quark mass, //c, 3 - This time, there is a visible discrepancy between rnD{|J^c)\^lc=^J.c ,3 
the corresponding linear extrapolation. This difference, however, is less than lOMeV, i.e. sig¬ 
nificantly smaller than the estimated magnitude of lattice discretization errors (cf. section [5]) . 
Moreover, the obtained valence charm quark mass /ic,phys) which is rather close to is essen¬ 
tially identical to a valence charm quark mass one would obtain, when using instead of 

/^c,l, /ic,2- 

Note that there is a mismatch of around 10% ... 20% of the strange and charm sea quark masses 
of our ensembles (represented by = 0.150 and /j-s = 0.190; cf. section [2TT) and the corre¬ 
sponding physical values: ~ 0.59 GeV^ and mo,sea ~ 2.1 GeV. Since strange 

and charm sea quarks are expected to have a rather small effect on hadron masses this slight 
mismatch should have negligible influence on the resulting meson spectra. 

Now that physical strange and charm valence quark masses are known, the linear extrapolation 
procedure is reversed; each D meson, Dg meson and charmonium mass m is computed at two 
or three pairs of valence quark masses, {fJ-g^i, fJ-c, 2 ), IJ-c,i) and {fJ-g, 2 , IJ-c, 2 ), and the linear 

extrapolation is then performed to (/is,physj A*c,phys) according to 


^phys 


m{nc,2) + (m{^c,i) 



/^c,phys /^c,2 
A^c,l Me,2 


( 22 ) 
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for D mesons and charmonium and according to 


^phys 


= ^J^c,2) - 'mifis,2, Aic,2) ) ^ 

/^s,l /^s ,2 

\ ( ( \ ( \\ A^c,phys ~ Me,2 

+ U^(Ms,2,Mc,i) - ^(M5,2,Mc,2) -- 

' Mc,l Me,2 


(23) 


for Dg mesons. Of course, it would have been possible to perform computations of these meson 
spectra using directly strange and charm valence quark masses (Ais,phys) Mc,phys)- We consider 
this linear extrapolation procedure, eqs. (1221) and (1231) . however, superior, because it allows a 
lot more flexibility during the final analysis. For example, one can easily change the value of 
the lattice spacing to investigate a possible source of systematic erroiH without the need for 
redoing propagator computations and contractions. Alternatively, one can even determine the 
lattice spacing in physical units by matching the obtained meson spectra with corresponding 
experimental data. Another application would be to crudely estimate systematic errors due 
to isospin breaking and electromagnetic effects using e.g. 2m|,± — and m£)± instead of 
~ "^-^0 and when determining the strange and charm valence quark masses lJ,s,phys 

and ^c,phys- We plan to investigate such issues and use this flexibility for an evolved error 
analysis in an upcoming publication, when we have meson masses available for several values of 
the lattice spacing. 

In Figure [2] (bottom) we show examples of linear extrapolations of various Dg meson masses 
((±,=f) twisted mass sign combinations, representations Ai,Ti,r 2 and V = ±) both in 
fig and in //c to physical strange and charm valence quark masses (ensemble A80.24). Again 
we compare with computations performed with a third strange and charm valence quark mass 
(blue points). We find excellent agreement demonstrating once more the validity of the linear 
extrapolations. Similar consistent results have been obtained for D meson and for charmonium 
masses and for both (±,q=) and (±,±) twisted mass sign combinations. 


4.4 Determination of statistical errors 

Statistical errors for D and Dg meson and charmonium masses (extrapolated to physical strange 
and charm valence quark masses as explained in the previous subsection) are determined on each 
ensemble via an evolved jackknife analysis starting at the level of the correlation matrices (1141) . 
To exclude statistical correlations between gauge link configurations, which are close in Monte 
Carlo simulation time, we performed a suitable binning of these configurations. 


5 Results 

We compute D meson, Dg meson and charmonium masses on each of the three ensembles listed 
in Tabledl These ensembles correspond to a single lattice spacing a ~ 0.0885 fm, but differ in the 
unphysically heavy u/d quark mass (m,r ~ 276MeV, 315 MeV, 443MeV). Each meson mass is 

®There appear to be unresolved inconsistencies between different collaborations regarding scale setting and 
standard non-perturbative scales like ro [B3]. Even within the ETM Collaboration there exist two values of the 
lattice spacing for the ensembles we are using: a = 0.0885(36) fm obtained from the pion decay constant [56] and 
a = 0.0920(21) fm obtained from the nucleon mass |64 |. 
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extrapolated linearly in to the physical value of the u/d quark mass (ttItt = rn^o = 135 MeV). 
Strange and charm valence quark masses correspond to their physical values (cf. section . 

Since at the moment computations are only available for a single lattice spacing, we are not 
in a position to perform a continuum extrapolation. Lattice discretization errors are, however, 
expected to be small, because the lattice spacing we use, a « 0.0885 fm, is rather fine and meson 
masses are automatically 0{a) improved due to our specific quark discretization, Wilson twisted 
mass lattice QCD at maximal twist |54] . In this formulation each meson mass can be computed 
in two slightly different ways, either using (±,q:) or (±,zb) twisted mass sign combinations (cf. 
section [3.2.ip . The corresponding two results differ by lattice discretization errors. We use these 
differences to crudely estimate the magnitude of the lattice discretization errors associated with 
our resulting meson masses. 

The continuum limit will be studied in a future publication. There we also plan to include a 
more detailed discussion and analysis of systematic errors, including e.g. isospin breaking and 
and electromagnetic effects. 

5.1 D and Ds mesons 

Since the D and the Dg meson spectrum are qualitatively very similar, these two sectors are 
discussed in parallel in the following. The charmonium sector is presented separately in sec¬ 
tion [5]2j 

5.1.1 Ai representation (spin J = 0) 

We compute the masses of the two lowest states in the Ai sector (both for D and for Dg mesons), 
one with V = —, the other with V = +, for each of the three ensembles (cf. Tabled]) and both 
twisted mass sign combinations (±,t) and (±,±). In this process we solve the generalized 
eigenvalue problem (|16l) using 2x2 correlation matrices with creation operator indices 1,3 (cf. 
Table [2|), i.e. T = 75,1. Note again that parity is not an exact symmetry in Wilson twisted 
mass lattice QCD at finite lattice spacing. The positive and negative parity states form a single 
sector, where the ground state has V = — (the D/Dg meson) and the first excited state has 

V = + (the Dq/D*q meson). One has to determine the masses of both states at the same time 
from a single correlation matrix (cf. the detailed discussion in section l4^ . 

Even though we have implemented eight different creation operators for the Ai representation 
(cf. Table [2|), the above mentioned rather small 2x2 correlation matrices turn out to be an 
optimal choice: statistical errors are quite small and long and stable plateaus are obtained both 
for effective masses and for operator contents (cf. e.g. Figure d])- 

V = —: D and Dg 

As already mentioned above and as expected the ground states in the Ai sectors have negative 
parity, i.e. correspond to the D meson and the Dg meson. In Figure [3] the computed D meson 
masses (left) and Dg meson masses (right) are shown for all three ensembles (i.e. 
tUtt ~ 276MeV, 315MeV, 443MeV) and for both twisted mass sign combinations, (±,t) (red 
points) and (±,±) (green points). The straight red and green lines are linear extrapolations in 
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to physical u/d quark masses corresponding to = m^o = 135 MeV. The results of these 
extrapolations are shown in magenta. The blue points are experimental results for m^jo, 
and rriDs (note that there are two experimental results in the left plot, because the neutral 
and the charged D meson differ in mass, m£,± — m£,o = 5 MeV)@ 
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Figure 3: Ai representation (spin J = 0), V = —. (left): D meson, (right): Dg meson. 


For each ensemble the valence charm quark mass is chosen such that the (±,t) version of 
the D meson mass is identical to the experimental result = 1865 GeV, as discussed in 

section 01 Therefore, these (±,t) D meson masses should not be considered as predictions 
and, consequently, a linear extrapolation of these masses (which would trivially be a constant) 
is neither needed nor shown. The statistical errors of these masses enter an evolved jackknife 
procedure and are, therefore, not only considered in the errors of the valence quark masses 
/ic,phys; but also in the errors of all D meson, Dg meson and charmonium masses computed in 
this work. 

The differences between the lattice QCD results obtained with (±,t) and with (±,±) twisted 
mass sign combinations are for both the D and the Dg meson around 50 MeV. These differences, 
which will vanish in the continuum limit, are a crude estimate of the magnitude of lattice 
discretization errors associated with our current results obtained at a single lattice spacing, i.e. 
relative errors of around 2.5%. Note that for meson masses obtained with (±,t) twisted mass 
sign combinations discretization errors are expected to be significantly smaller ISIIES], i.e. the 
mentioned 50 MeV are most likely a rather conservative estimate. This is consistent e.g. with our 
(±,7p) lattice result for the Dg meson, which differs by less than 10 MeV from the experimental 
result. 

The linear increase of the Dg meson mass for decreasing u/d quark mass is an expected con¬ 
sequence of our procedure for setting the charm valence quark mass /ic,phys- We choose ^c,phys 
for each ensemble such that the lattice result for mo agrees with the experimental result 

®The majority of plots shown in this section, including the two plots in Figure [S] follows this style. The right 
panel of each plot is a zoomed version of the left panel with respect to the vertical axis, but otherwise identical. 
In the left panels the scale of the vertical axis, which represents the meson mass, always ranges from 0 to 4.6 GeV. 
Hence from the left panels one can conveniently read off the relative errors and precision of our results. They also 
allow to directly compare meson masses from different plots. The right panels are strongly zoomed (individually 
for each meson) and, hence, show more clearly details regarding the data quality, the absolute size of the errors 
and the dependence of the meson mass on the u/d quark mass. 
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m^jo = 1865 MeV, i.e. independently of the u/d quark mass (cf. section 031) ■ Clearly an in¬ 
creasing u/d quark mass leads to a decreasing /Uc,phys slightly lighter than the physical charm 
quark mass. This in turn yields the observed u/d quark mass dependence of rrif)^. 

The experimentally observed splitting m£,± — 111^0 = 5 MeV indicates the magnitude of elec¬ 
tromagnetic and isospin breaking effects. Since the currently estimated discretization errors of 
<50 MeV are much larger, we will at the moment not consider systematic errors due to the 
neglect of electromagnetism and isospin breaking. We plan to address such errors in a future 
publication, where we will study the continuum limit. 


V = +: Dl and 

The first excitations in the Ai sectors have positive parity, i.e. should correspond to the Dq 
meson and the D*q meson (cf. Figure 0]). While the masses of these states have been extracted 
from a 2 X 2 correlation matrix with creation operators T = 75 and T = 1 (indices 1 and 
3 in Table [2]) , we have studied their structure by considering also larger correlation matrices 
with additional V = + creation operators (indices 4, 7 and 8 in Table [2|). Interestingly, the 
operator contents of both the Dq meson and the D*q meson are mixtures of T = 1,70 and 
T = 7 jnj, 7 o 7 jnj of roughly the same magnitude. This indicates that the quarks inside these 
mesons form superpositions of S and P waves and not predominantly P waves as suggested by 
many quark models. 


Dq meson 



in GeV^ 


in GeV^ 



meson 



; in GeV' 


Figure 4: Ai representation (spin J = Q),V = +. (left): Dq 


meson. 


(right): D 


sO 


meson. 


In both cases we observe decreasing meson masses for decreasing u/d quark mass, a behavior not 
present for the previously discussed V = — parity partners. An explanation could be that the 
extracted V = + states contain rather light two-meson contributions with the same quantum 
numbers, D + tt and D + K, respectively. In Figure 0] this is illustrated by the gray curves, which 
correspond to the estimated masses of the two-meson states, mn + and mo + mKinT-w) 0 - 
For the Dq meson this is also supported by the rather large x^/doi of the extrapolation to 
physical u/d quark mass (cf. Figure 0] [left]). Surprisingly such a behavior has not been observed 
in m, where a rather similar lattice QCD setup, in particular a similar set of creation operators, 
has been used. 

^Here and in the following we have used mx(m- 7 r) = -\/(0.477 GeV^ -|- m|)/2, mn = 1.865 GeV and m_D» = 
2.007 GeV for these estimates. 
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Since it is unexpectedly light, in the literature the D*q meson is frequently discussed as a mesonic 
molecule or tetraquark candidate (cf. e.g. [65l[66]). In principle it is possible to investigate the 
structure of these states using lattice methods, but this will require the implementation of 
additional four-quark creation operators of mesonic molecule, of diquark-antidiquark and/or of 
two-meson type. We are in the process of developing techniques for computing corresponding 
correlation matrices [38l|39lll0]. Recent lattice papers using four-quark operators and focusing 
specihcally on Dq and are [28l [29l |30l (3^ 1^ . 

This time lattice discretization errors indicated by the differences between results obtained 
with (±,t) and (±,±) twisted mass sign combinations are somewhat larger, of the order of 
60MeV ... 100 MeV (i.e. relative errors ~ 3% ... 5%). Within this crude conservative estimate 
there is consistency with experimental results. It is reassuring that the (±,t) result for m£)*, 
which is expected to have signihcantly smaller discretization errors than the (±,±) result, is 
within Ifj of the experimental result for ( 2400)0 without taking discretization errors into 
account. The (±,t) Dsq result is around 80MeV larger than its experimental counterpart. 
This is similar to what has been found in quark models (cf. e.g. [5]) and other lattice QCD 
computations using exclusively quark-antiquark meson creation operators (cf. e.g. [E]). It could 
be an indication that the D*q meson is not predominantly a quark-antiquark state, but possibly 
a mesonic DK molecule (e.g. supported by m) or a diquark-antidiquark pair. 

5.1.2 Ti representation (spin J = 1) 

V = -: D* and D* 

We proceed as in section 15.1.11 where we have discussed the Ai representation. This time we 
have solved generalized eigenvalue problems (|16l) using 4x4 correlation matrices with creation 
operator indices 1,2,3,4 (cf. Table [2]) including the commonly used V = — operators T G 
{ 71 , 7071 }. As expected the ground states in the Ti sectors have negative parity, i.e. correspond 
to the D* meson and the D* meson (cf. Figured]). 


D meson Dg meson 




Figure 5: Ti representation (spin J = 1), R = —. (left): D* meson, (right): D* meson. 

While at physically light u/d quark masses the decay of D* to H -|- vr is possible, it is excluded 
in our computations with pion masses m7r>276MeV. This is consistent with the rather mild 
dependence of tud* (cf. Figure [5] [left]). 
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The differences between the lattice QCD results obtained with (±, and with (±, ±) twisted 
mass sign combinations are around 10 MeV ... 20 MeV, i.e. lattice discretization errors seem to be 
smaller than in the Ai representation. There is also a discrepancy of around 20 MeV ... 40 MeV 
to the corresponding experimental results [T]. In this context it is interesting to note that the 
mass differences ra£)* — rajy and mu* — mus are rather strongly dependent on scale setting, i.e. 
the lattice spacing in physical units. For example, when using a = 0.0920(21) fm (obtained from 
the nucleon mass [64]) instead of a = 0.0885(36) fm (obtained from the pion decay constant 
|56ji. the above mentioned discrepancy of 20 MeV ... 40 MeV is reduced to 15 MeV ... 20 MeV. 
We consider this as an indication that the larger lattice spacing obtained from the nucleon mass 
might be better suited to determine hadron masses with small discretization errors than the 
lattice spacing obtained from the pion decay constant (which is more commonly used within 
the ETM Collaboration). Similar observations have been reported in spectrum computations 
of B and Bg mesons [591 EZ] and b baryons [68]. As mentioned before we plan to investigate 
such scale setting issues and the continuum limit in detail in a future publication, when we have 
meson masses available for several values of the lattice spacing. 


V = +: T»i(2430), T>i(2420) and T»si(2460), T>si(2536) 

To determine masses of positive parity states we use larger 8 x 8 correlation matrices with 
creation operator indices 1, 2,3,4, 7, 8 ,9,10 (cf. Table|2|), i.e. the previously used V = — operators 
T G { 71 , 7071 } and all available V = + operators with L = 0 and L = 1. The first and second 
excitations in the Ti sectors have positive parity, i.e. correspond to i7i(2430), T)i(2420) and 
77,i(2460), T>,i(2536). 

The masses of the two states in each sector are rather close. Therefore, it is not obvious, 
how to correctly assign the obtained lattice QCD results to i7i(2430) and i7i(2420) and to 
Dsi(24:60) and i7si(2536), respectively. Such an assignment is not only important for a complete 
and precise computation of the D and Dg meson spectrum, but also in the context of specific 
decays, in particular B^*'> —>• D** + l + where D** denotes the four positive parity D mesons 
with J = 0,1,2, which include Di(2430) and i7i(2420). To understand these decays is e.g. 
essential for a precise determination of the CKM matrix element Vcb- However, there is a long¬ 
standing conflict between theory and experiment regarding the corresponding branching ratios 
m- While these decays have been studied with lattice QCD in the static limit some time ago 
|69j , recently computations with b and c quarks of finite mass have been started |42l 143] . The 
latter computations are, however, restricted to D** with J = 0, 2, mainly because the separation 
of the two J = 1 states is rather difficult. In the following we demonstrate, how to distinguish 
those two states and correctly assign their masses to Di(2430) and Di(2420) (and similarly to 
Dsi(2460) and Dsi{2536) in the Dg sector). 

At first it is important to note that even though Di(2430) and Di(2420) have similar masses, 
their structure is quite different. Due to the heavy charm valence quark D mesons are expected 
to be qualitatively similar to static-light mesons. Since the spin of a static quark is irrelevant, 
it is appropriate to label static light mesons by the half-integer total angular momentum j of 
their light degrees of freedom, i.e. the light quarks and gluons. One of the two states Di(2430) 
and Di(2420) has j ~ 1/2, while the other has j ~ 3/2 (this expectation has been confirmed by 
model calculations, e.g. M)- For a detailed discussion cf. e.g. [591 [67]. 

The experimental results can be classified according to j ~ 1/2 and j ~ 3/2 as follows: 
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• D mesons: 

Both L)i(2430) and -Di(2420) can decay to D* + tt. Di( 2430) has a rather large width 
(r = 384t|ioMeV), whereas L)i(2420) is comparably stable (F = 27.4 ± 2.5 MeV). This 
difference in the widths suggests the assignment j « 1/2 to Hi (2430) and j « 3/2 to 
Hi (2420): while Hi (2430) can then readily decay via an S wave, Hi (2420) is protected by 
angular momentum j ~ 3/2, which only allows a less likely H wave decay. 

• Ds mesons: 

Only Hsi(2536) can decay to D* + K. Hsi(2460) is too light for such a decay, since 
™Dsi( 2460 ) < Both states, however, have rather small widths, F < 3.5 MeV 

and F = 0.92 ± 0.05 MeV, respectively. Consequently, the heavier state Hsi(2536) must 
be protected by angular momentum, i.e. have j « 3/2, while the other state Hsi(2460) 
corresponds to the remaining j ~ 1/2. 


To decide, which of the two lattice QCD results for masses of = !■*■ H meson states cor¬ 
responds to j ~ 1/2 and which to j ~ 3/2, we study the eigenvectors obtained by solving the 
generalized eigenvalue problem (jlOp . We use linear combinations of Ti meson creation operators 
from Table [21 which excite not only states with definite J, but also with definite light total 
angular momentum j. These 16 linear combinations are collected in Tableland sorted into five 
classes Ci,..., Cs, one with V = — and four with V = +, where the latter differ inj = 1/2,3/2 
and L = 0,1,2. For each extracted mass n = 0,1,2 (0: V = — ground state; 1,2: V = + 
excitations) and for each class k = 1,..., 5 we sum over the squared components of the vectors 
(cf. eq. (HH)), 


zt\t) = 


E 

lilOiGCfe} 




(24) 


(n) 

These quantities Z/ ^ indicate the parity of the state n and in case oi V = + its light total 
angular momentum j and orbital angular momentum L. 

obtained by solving 16 x 16 generalized eigenvalue problems are plotted in Figure |6| as 
functions of the temporal separation t for H mesons and the A30.32 ensemble. The upper 
plot confirms that the ground state has V = —, i.e. corresponds to D*. The lower left plot 
corresponding to the first excitation shows a strong dominance of {V = +,j = 1/2, L = 0) and 
{V = +,j = 1/2, L = 1). Clearly, this state should be interpreted as j « 1/2, i.e. as Hi(2430), 
where L = 0 has a somewhat larger contribution than L = 1. The second excitation (lower right 
plot) is almost exclusively j gs 3/2, i.e. corresponds to Hi (2420), with L = 1 (this is the reason, 
why we discard L = 2 meson creation operators in the final determination of = 1+ H meson 
masses; cf. also the first paragraph of this subsection). 

In this context we also refer to the lattice QCD study [29], where the identification of j ~ 1/2 
and j ~ 3/2 has been done applying a different strategy, namely including and excluding meson 
creation operators composed of four quarks. One of the two = 1+ states is essentially 
unaffected and, therefore, interpreted as the stable Hi (2420) meson with j ~ 3/2, while the 
other state exhibits a certain sensitivity with respect to the inclusion/exclusion of four quark 
creation operators and, hence, is interpreted as the less stable Hi (2430) meson with j ~ 1/2. 

We have carried out a similar analysis for = 1^ Dg mesons and obtained qualitatively 
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class 

continuum 

twisted mass lattice QCD 

V 

j 

L 

r(n),pb 

tb, (±,t) 

tb, (±,±) 

Cl 

- 

1/2 

0 

7i 

7o7i 

±*75 X 
pb 

pb 

±*75 X 

1 

75 (n x^)j+ iyoiij 
7o75(n X 7) ,- -F n,' 

pb 

±*75 X 

±*75 X 

pb 

3/2 

1 

75(n X 7)j - 2z7onj 
7o75(n X 7),- - 2nj 

±*75 X 

pb 

pb 

±*75 X 

2 

7i(2n^ -n^ -ng) 
707i(2n? -n2 -n2) 

±*75 X 

pb 

pb 

±*75 X 

C2 

+ 

1/2 

0 

7571 

707571 

±*75 X 

pb 

pb 

±*75 X 

C3 

+ 

1/2 

1 

(n X 7 )j - i7o75nj 

70(n X 7)j - 

±*75 X 

pb 

pb 

±*75 X 

C4 

+ 

3/2 

1 

(n X 7)j + 2i7o75nj 

70 (n X 7)j -b2i75nj 

pb 

±*75 X 

±*75 X 
pb 

Cs 

+ 

3/2 

2 

7571 (2n^ -n^ -n^) 
707571 (2nf - - n|) 

±*75 X 

pb 

pb 

±*75 X 


Table 4: = I"*" meson creation operators with definite j and L. 


identical results: the lighter of the two extracted states has j « 1/2 (i.e. corresponds to the 
Zlsi(2460) meson), while the heavier has j ~ 3/2 (i.e. corresponds to the Zlsi(2536) meson). 

The J'^ = 1+ D and Ds meson masses are shown in separate plots for j ~ 1/2 and j ~ 3/2 in 
Figure [71 

The lattice result for the mass of Hi (2430) (j ~ 1/2; top left plot) is in perfect agreement with 
the corresponding experimental result. Nevertheless it should be treated with caution. Since 
the Hi (2430) has a large width and is expected to be rather unstable (it can decay in D* + vr), 
a solid and more rigorous result will require a proper resonance treatment as e.g. pioneered in 

m- 

The corresponding Dg mass with j ~ 1/2 (Hsi(2460); top right plot) is about 70MeV larger than 
its experimental counterpart. As before for H*q, this is similar to what has been found using 
quark models, e.g. [5]. It could be an indication that the Hsi(2460) meson is not predominantly 
a quark-antiquark state, but possibly a mesonic D*K molecule or a diquark-antidiquark pair. 

The lattice QCD results for the masses of the j « 3/2 states. Hi (2420) and Hsi(2536) are also 
around 50 MeV ... 100 MeV larger than the corresponding experimental results. This mismatch 
is somewhat surprising, since both states are rather stable and quark model calculations assum¬ 
ing a straightforward quark antiquark structure are able to reproduce the experimental values 
quite accurately. We presume that this discrepancy of around 2% ... 4% is due to discretization 
errors. For example when using the lattice spacing a = 0.0920(21) fm obtained from the nucleon 
mass |64] . the discrepancy is only half as large. It will be interesting to see, whether there will 
be agreement with experimental results, after performing a continuum extrapolation. 
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operator content of the ground state (n = 0) 



operator content of the first excitation (n = 1) 


operator content of the second excitation (n = 2) 
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- 

C 4 : P = +, j = 3/2, L = 1 . 



C 5 : P = +, j = 3/2, L = 2 . 






. 



Figure 6: Operator content for = 1 and D mesons (A30.32 ensemble). 


5.1.3 E and T 2 representations (spin J = 2) 

Effective mass plateaus obtained by solving generalized eigenvalue problems (jl6[) are rather short 
for the E and the T 2 representation, i.e. the determination of meson masses is more subtle than 
for the Ai and the Ti representation before. 

There seem to be stronger statistical fluctuations for the E representation than for the T 2 
representation. Therefore, for the E representation we only determine the ground state, which 
has V = +. We use 4x4 correlation matrices with creation operator indices 5 to 8 (cf. Table [2]), 
i.e. operators with angular momentum L = 1. 

For the T 2 representation it is possible to extract additionally two V = — states. The corre¬ 
sponding correlation matrices contain creation operators with indices 1 to 4 (operators with 
L = 1; cf. Table [2]) and for the Dg sector also operators with indices 5 and 6 (operators with 
L = 2; cf. Table [2]). 

V = +: £>^(2460) and T>*2(2573) 

The ground states in the E and the T 2 representations have positive parity and are of similar 
mass. This strongly indicates that these are J = 2 states, since total angular momentum J = 2 
is part of both E and T 2 , in contrast to e.g. J = 3. Consequently, we interpret them as 112(2460) 
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Figure 7: Ti representation (spin J = 1),V = +. (top left): F)i(2430) meson {j « 1/2). (top 
right): Zlsi(2460) meson (j ss 1/2). (bottom left): Zli(2420) meson {j « 3/2). (bottom 
right): Zlsi(2536) meson {j « 3/2). 


and 71*2(2573). 

Within statistical errors the meson masses are essentially independent of (cf. Figure [8]). This 
indicates rather stable states and is expected, because total angular momentum J = 2 allows 
decays 112(2460) —)• H + vr and 11*2(2573) D + K, respectively, only via D waves. Such decays 
are strongly suppressed compared to S wave decays, which are possible e.g. for corresponding 
= O"'' and (J^ = ~ 1/2) states (cf. the more detailed discussion in section 15.1.21) . 

Discretization errors indicated by the differences between results obtained with (±, ip) and (±, ±) 
twisted mass sign combinations as well as from the E and the T 2 representations are together 
with statistical errors of the order of 100 MeV (i.e. relative errors ~ 4%). The extrapolation 
to physically light u/d quark mass yields meson masses, which are also around 100 MeV larger 
than the corresponding experimental results. As discussed for the (J = 1, j = 3/2) states in the 
previous subsection, discretization errors might be the reason for this discrepancy. 


V = -■. 11(2750) 

Resulting meson masses from the T2 representation are shown in Figure [9l Crude results with 
large statistical errors from the E representation (not shown in Figure [9]) are in agreement and, 
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Figure 8: E and T 2 representations (spin J = 2), V = +. (left): 1)2(2460) meson, (right): 
F»*2(2573) meson. 


hence, suggest J = 2. There seems to be only little dependence on indicating rather stable 
states. Note, however, that these negative parity states are less strongly protected by angular 
momentum than their parity partners 1)2(2460) and Z)*2(2573), since P wave decays to D* + vr 
and D* + K are possible. 
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Figure 9: T 2 representation (spin J = 2), P = —. (left): D(2750) meson, (right): Dg sector. 

The experimentally observed 7)(2750) is usually interpreted as a = 2~ state. We find 
agreement with this experimental result [T] within statistical errors as well as with recent lattice 
QCD computations of the mass of the D meson with = 2~ (at unphysically heavy u/d 
quark mass corresponding to m-,^ ps 266 MeV and m-,^ ps 391 MeV, respectively) [29l[23]. The 
corresponding = 2~ Dg state has experimentally not yet been measured or clearly identified. 
Therefore, our lattice result can be considered as theoretical prediction. We find meson masses 
consistent with the lattice QCD computation [21] and quark model predictions from [5|. 


5.2 Charmonium 

We compute masses of charmonium states using both twisted mass sign combinations (±, ^p) 
and (±,±) neglecting disconnected contributions to the correlation matrix elements (I14|] . It 


25 





























is expected that this neglect of disconnected diagrams introduces only tiny systematic errors, 
which are much smaller than our combined statistical uncertainties and lattice discretization 
errors, because of the OZI suppression of these diagrams. This expectation is quantitatively 
supported e.g. by a quenched lattice QCD computation, where the effect of disconnected dia¬ 
grams on the charmonium hyperfine splitting is found to be ~ 1... 4 MeV m- A corresponding 
consistent perturbative estimate is a shift of charmonium masses by ~ 2.4 MeV (cf. |71l 1721 [73] 
and references therein). 

Note that twisted mass quantum numbers are different, when using (±,t) and (±,±) twisted 
mass sign combinations. For (±, ±) charge conjugation C is a quantum number (but not parity 
V), while for (±, only the product P o C is a quantum number (cf. also section [4T1 where we 
list the sizes of our correlation matrices, or, alternatively, Tabled columns “C” and “PC”, re¬ 
spectively). All charmonium analyses are based on the maximal sets of available meson creation 
operators, i.e. 2 x 2, 4 x 4, 6 x 6 or 10 x 10 correlation matrices as listed in section El 

5.2.1 Ai representation (spin J = 0) 

V = -•. r]c{lS) and r]c{2S) 

In Figure [ini (left) = 0 charmonium masses are shown. 


r|(,(1S) and r|;,(2S) mesons 
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Figure 10: Ai representation (spin J = 0). (left): VC = — h, ??c(l*S') and r]c{2S) mesons, 
(right): VC = -h-h, Xco(lP) meson. 

The ground state and the first excitation can clearly be identihed with the experimentally known 
r]c{lS) meson and r]c{2S) meson. Even though we use a heavy c quark and a heavy c antiquark, 
discretization errors indicated by the differences between results obtained with (±, =p) and (±, ±) 
twisted mass sign combinations are rather small, around 70 MeV, i.e. ~ 2.5%. Moreover, within 
this uncertainty there is perfect agreement with experimental results [T]. The fact that our lattice 
QCD result for the mass of r]c{lS) at finite lattice spacing is slightly below the experimental 
result is also in agreement with a 2-flavor twisted mass lattice QCD study including the r]c{lS) 
meson m- 

We are also able to extract a crude signal for a third state, for which there is currently no 
clearly identified experimental counterpart. Such a state has, however, been observed in previous 
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lattice QCD computations, e.g. in [20l[29] (at unphysically heavy u/d quark mass corresponding 
to TTijr « 266 MeV and ps 391 MeV, respectively), around 400 MeV... 500 MeV above the 
mass of the r/c(2S) meson. Here we observe an unexpectedly large difference of around 200 MeV 
between (±,=p) and (±,±) results. The reason for this could be the rather bad plateau quality 
of the corresponding effective masses (short plateaus formed by only three points, before the 
signal is lost in noise), which might be a sign of contamination by even higher excitations. This 
is also supported by the fact that our result is around 600 MeV... 800 MeV heavier than the 
mass of r]c{2S), i.e. heavier than the above mentioned lattice QCD results from [201129] . Another 
reason could be mixing with a lighter state of the VC = ++ sector: for the (±, ±) result we 
observe a rather large contribution of around 30% of the operator r(n) = 'jQjjiij corresponding 
to quantum numbers VC = ++ (index 8 in Table [2]), while for (±,t) such a mixing is excluded 
by the symmetry V oC. 

Even though we use a meson creation operator with quantum numbers VC =-(index 5 in 

Table [2]), we do not obtain a clear signal, i.e. a trustworthy effective mass plateau, to extract a 
corresponding meson mass. 


V = +: XcoilP) 

In Figure do] (right) the resulting masses for Xco{^P) (quantum numbers = O^"*") are shown. 
Discretization errors indicated by the differences between results obtained with (±, ^) and (±, ±) 
twisted mass sign combinations are rather small, around 50MeV, i.e. ~ 1.5%. The u/d quark 
mass extrapolated (±, t) result, for which one expects smaller discretization errors, is in perfect 
agreement with its experimental counterpart [I]. 

Similarly to the VC =-sector we are not able to reliably determine an exotic VC = 4— state. 

5.2.2 Ti representation (spin J = 1 and J = 3) 

V = -: J/T(15), T(25) and T(3770) 

For VC =-we are able to extract the masses of four states, which are shown in in Figure dH 

(top). 

The two lowest states can be identified as the J/T(1S') and the 'I'(25) meson. The lattice 
QCD results for both twisted mass discretizations (±,t) (±,±) and the corresponding 

experimental results [I] differ by less than 50 MeV, i.e. are in agreement within the expected 
discretization errors. 

The second and third excitation are very close in mass, around 3800 MeV... 3900 MeV. One of 
these two states should correspond to the T(3770) meson. The other state seems to correspond 
to a J = 3 state in the continuum (the continuum spins of Ti are J = 1,3,4,...), because a 
state with the same mass is observed in the T 2 representation (continuum spins J = 2, 3,4,...), 
but not in the E representation (continuum spins J = 2,4,...; cf. section r5.2.3p . A = 3 
state in this energy region has also been observed in other lattice QCD studies (cf. e.g. [201129]) 
and has also been predicted using Dyson-Schwinger and Bethe-Salpeter equations [7]. 

Even though we use two creation operators with exotic quantum numbers VC = —h (indices 6 
and 12 in Table [2]), we do not obtain a clear signal, to reliably extract a corresponding mass. 
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Figure 11: Ti representation (spin J = 1 and J = 3). (top): VC = -, J/'I'(15), ^^(25) and 

'I'(3770) mesons and a state with J = 3. (bottom left): VC = ++, Xci(l-P) meson, (bottom 
right): VC = H—, hc{lP) meson. 


V = +: Xci(l-P) and hc{lP) 

In Figure [11] (bottom) the resulting masses for Xci(l-P) (quantum numbers = 1"’"*') and 
hc{lP) (quantum numbers = I"' ) are shown. The lattice QCD results for both twisted 
mass discretizations (±,t) and (±,±) and the corresponding experimental results [T] differ by 
less than 40MeV, i.e. are in agreement within the estimated discretization errors. 


5.2.3 E and T 2 representations (spin J = 2 and J = 3) 

V = +• Xc 2 (l-P) and Xc2(2P) 

For V = -\- we are able to extract two states with C = + both for the E and the T 2 represen¬ 
tation (cf. Figure fT^ [top]). Since the masses of the two ground states as well as the masses of 
the two excitations agree within errors, we interpret them as J = 2 states, i.e. as the Xc 2 i^P) 
meson and the Xc2(2F’) meson. Within statistical errors there is essentially no dependence on 
m^. Discretization errors indicated by the differences between results obtained with (±,t) and 
(±,±) twisted mass sign combinations as well as from the E and the T 2 representations are for 
the Xc 2(1F) meson together with statistical errors of the order of 50MeV (i.e. relative errors 
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~ 1.5%). Within this estimated uncertainty there is excellent agreement with the corresponding 
experimental result. For the Xc 2 {‘^P) meson the uncertainty is roughly twice as large, i.e. around 
100MeV. Our lattice QCD result is also slightly larger (around 150MeV i.e. l.ba larger) than 
the experimentally measured mass, which could be an indication that there is a small contami¬ 
nation of higher excitations (cf. also the lattice QCD study [20] . where a similar trend has been 
observed). 


Xc2(1 P) and Xc 2(2P) mesons 
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Figure 12: E and T 2 representations (spin J = 2 and J = 3). (top): VC = Xc 2 i^P) and 

Xc 2(2-P) mesons, (bottom left): VC = —+, r]c 2 meson, (bottom right): VC = -, one state 

with J = 2 and (in the T 2 representation) another state with J = 3. 


V = -: r/c2 

For V = — we are able to determine one state with C = + both for the E and the T 2 represen¬ 
tation, which should correspond to the 7]c2 meson (cf. Figure [12] [bottom left]). Again combined 
discretization and statistical errors are around 50 MeV. Currently there is no established exper¬ 
imental result available, to which we can confront our result. We find, however, agreement with 
theoretical predictions from other lattice QCD computations [201129] and from a calculation 
using Dyson-Schwinger and Bethe-Salpeter equations [7]. 

For V = — and C = — in the mass region 3800 MeV... 3900 MeV we find only a single state for 
the E representation (continuum spins J = 2,4,...), but two states for the T 2 representation 
(continuum spins J = 2,3,4,...) (cf. Figure ITT] [bottom right]). This suggests to interpret the 
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state, which is present both in the E and the T 2 representation, as a J = 2 state in the continuum 
and the additionai state in the T 2 representation as a J = 3 state in the continuum. This is 
supported by our resuits for the Ti representation (continuum spins J = 1,3,4,...), where we 
have extracted a state in the same mass region (cf. section [5.2.21 in particuiar Figure [H] [top]). 
Again we can oniy compare to other theoreticai predictions, e.g. from [221ESE!, which are in 
agreement with our resuits. 


6 Summary and conclusions 

We have computed masses of iow iying D meson, Dg meson and charmonium states with totai 
anguiar momentum J = 0,1,2,3, parity V = —, + and charge conjugation C = —, + using 
Wiison twisted mass lattice QCD. We have used gauge link ensembles generated by the European 
Twisted Mass Collaboration with three different u/d quark masses corresponding to € 
{276MeV, 315 MeV, 443MeV}. After performing computations on these three ensembles, we 
have extrapolated the resulting meson masses to physically light u/d quark mass. 

Our computations are currently limited to a single lattice spacing, a ~ 0.0885 fm. Therefore, we 
are not able to study the continuum limit at the moment. In Wilson twisted mass lattice QCD it 
is, however, possible to compute meson masses using two different discretizations, either (±, =p) 
or (±,±) twisted mass sign combinations in the meson creation operators (cf. section [32]) • The 
differences between the resulting meson masses can be considered as crude estimates of lattice 
discretization errors. For the majority of mesons we have found differences of the order of 
50 MeV, which we take as an estimate of discretization errors. 

We expect these discretization errors to be the currently dominating source of systematic un¬ 
certainty. Further sources of error are listed in the following: 

• Finite spatial volume: 

Since for all three ensembles the spatial volume is rather large, i.e. rriT^L > 4, we expect that 
finite volume corrections are negligible compared to the above mentioned discretization 
errors of ~ 50 MeV. 

• Disconnected diagrams: 

When computing correlation matrices for charmonium states, we have omitted discon¬ 
nected contributions. These contributions are expected to be very small due to OZI sup¬ 
pression. The corresponding systematic errors for charmonium masses have been estimated 
to be less than 4 MeV (cf. section 15.21 and (TO] Ell ESI [73| ). 

• Electromagnetism and isospin breaking: 

We estimate the magnitude of electromagnetic corrections and effects due to isospin break¬ 
ing by comparing experimental results for masses of essentially stable charged and neutral 
D mesons (D, D*, Z)i(2420), D 2 ), yielding corresponding systematic errors < 5 MeV. This 
is consistent with [ZHEllEal, where such effects have been estimated to be ~ 2.6 MeV 
using a potential model. 

• Extrapolations in the up/down quark mass: 

The extrapolations to physically light u/d quark mass are linear in m^. To quantify a pos¬ 
sibly associated systematic error, one could compare different strategies of extrapolation. 
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e.g. also including a quadratic term in or using parameterizations of meson masses ob¬ 
tained from effective theories. We anticipate the corresponding uncertainty to be smaller 
than the currently estimated discretization errors of around 50MeV. We plan to study 
this issue in more detail in a future publication, where we will address the continuum limit. 

• Tuning of strange and charm sea quark mass: 

A probably tiny error might arise due to the slight deviation of strange and charm sea 
quark masses from their physical values (cf. section 14.3p . 

The final results of this work for D meson, Dg meson and charmonium masses are shown in 
Figure fOl and collected in Table EJ The numbers we quote have been obtained in the following 
way: 


• We use the u/d quark mass extrapolated results corresponding to (±, 4 :) twisted mass sign 
combinations, which are supposed to exhibit smaller discretization errors than their (±, ±) 
counterparts [ 6 II [62] . 

• At the moment systematic errors are expected to be strongly dominated by discretization 
errors (cf. the detailed discussion above). Therefore, we estimate all systematic errors to 
be around or less than 50MeV, the previously mentioned typical difference between meson 
masses computed with (±,t) and (±,±) twisted mass sign combinations. Since we use 
(±, 4 =) results, we consider this estimate to be rather conservative. 

• We assume independence of statistical and systematic errors. The total error is, hence, 
obtained by adding statistical errors and the 50MeV representing systematic errors in 
quadrature. 

• For J = 2 charmonium states we take the results from the E representation, where a mixing 
with and contamination by J = 3 states is excluded due to cubic rotational symmetry (in 
contrast to results from the T 2 representation). For D mesons and Dg mesons we take 
the results from the T 2 representation, which exhibit smaller statistical errors than the 
corresponding results from the E representation. 

• For the = 3 charmonium state we take the result from the Ti representation, which 
has smaller statistical errors than the corresponding result from the T 2 representation. 

In Figure [13] P = — states are shown in blue, while V = + states are shown in red. Statistical 
errors are represented by dark blue and dark red boxes, while the combined statistical and 
systematic errors are represented by light blue and light red boxes, respectively. The relative 
combined statistical and systematic errors of our results are in most cases between 2% and 3%. 
Within these errors our lattice QCD results agree with available experimental results [T], which 
are shown in gray. 

One of the next steps will be a computation of the same meson masses on several gauge link en¬ 
sembles at finer lattice spacings. This will enable us to perform a continuum extrapolation. Since 
the combined statistical and systematic errors are currently dominated by lattice discretization 
errors (crudely estimated by 50MeV), we expect that a continuum extrapolation will lead to 
significantly more precise results. 
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J'P 



D 

mesons 




Ds 

mesons 


lat. 

error 

exp. 

PDG 

lat. 

error 

exp. 

PDG 


mass 

sta./ 

com. 

mass 

name 

mass 

sta./com. 

mass 

name 

0- 

1865 



1870(0) 

1865(0) 


1962 

3 

50 

1968(0) 

Ds 

0+ 

2278 

20 

54 

2403(40) 

2318(29) 

Zlo*(2400)=^ 

2407 

12 

51 

2318(1) 

ll/o(2317) 

1- 

2051 

6 

50 

2010(0) 

2007(0) 

D*{2010)^ 

D*{2007)° 

2145 

4 

50 

2112(0) 

D* 

1+ 

2402 

23 

55 

2423(2) 

2421(1) 

Di{2A20)^ 

T»i(2420)° 

2521 

18 

53 

2460(1) 

Ai(2460) 

1+ 

2474 

35 

61 

2427(40) 

T>i(2430)° 

2608 

23 

55 

2535(0) 

ll^i(2536) 

2+ 

2549 

42 

65 

2464(2) 

2463(1) 

D2*(2460)=^ 

D|(2460)° 

2680 

18 

53 

2572(1) 

11 / 2 ( 2573 ) 

2“ 

2909 

66 

83 

2761(5) 

11(2750) 

2969 

33 

60 



2“ 

2996 

106 

118 



3001 

46 

68 





charmonium 


lattice 

error 

exp. 

PDG 


mass 

sta./com. 

mass 

name 

0-+ 

2884 

4 

50 

2984(1) 

Vc{lS) 

0-+ 

3535 

23 

55 

3639(1) 

Vc(2S) 

0-+ 

4372 

54 

74 



0++ 

3412 

6 

50 

3415(0) 

Xco(l-P) 

1“ 

3056 

5 

50 

3097(0) 


1— 

3662 

24 

56 

3686(0) 

i’i2S) 

1— 

3865 

23 

55 

3773(0) 

tp{3770) 

1++ 

3480 

9 

51 

3511(0) 

XciilP) 

1+- 

3494 

11 

51 

3525(0) 

hc{lP) 

2++ 

3536 

14 

52 

3556(0) 

Xc2{lP) 

2++ 

4066 

41 

64 

3927(3) 

Xc2i2P) 

2"+ 

3863 

21 

54 



2— 

3841 

12 

52 



3 — 

3865 

23 

55 




Table 5: Summary of lattice QCD results for D meson, Dg meson and charmonium masses (error 
sta.: statistical error; error com.: combined statistical and systematic error) and comparison to 
experimental results from [T]. 
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Figure 13: Summary of lattice QCD results for D meson, Dg meson and ctiarmonium masses 
{V = — states are shown in blue, V = + states in red; statistical errors are represented by dark 
blue/dark red boxes, combined statistical and systematic errors by light blue/light red boxes). 
For easy comparison experimental results from [T] are shown in black. 


Our strategy of computing meson masses using quark-antiquark creation operators allows to 
obtain solid and accurate results for states, which are mainly composed of a quark and an 
antiquark and which are quite stable. States, which might not fulfill these requirements, e.g. 
the rather unstable Dq, which readily decays into D + or I1 *q, which is frequently discussed 
as a tetraquark candidate, should finally be treated with more advanced lattice techniques. 
Further creation operators composed of four quarks (e.g. of mesonic molecule type, of diquark- 
antidiquark type and of two-meson type) have to be included in the correlation matrices. In case 
of unstable mesons, corresponding resonance parameters (mass, width) can then be extracted 
from the volume dependence of the spectrum of scattering states. We are currently in the process 
of developing and implementing such methods using a similar lattice setup [38l EU 00] . 

Another important aspect of this work is the separation and classification of the two = I"*" 
D meson states. Even though they have identical quantum numbers, their structure is quite 
different: one of them, F)i(2430), has j k, 1/2, while the other, F)i(2420), has j k, 3/2, where j 
denotes the spin and angular momentum of the light quark and gluons. Extracting those two 
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states unambiguously from a single = I"*" correlation matrix (both masses and eigenvector 
components, where the latter provide suitable linear combinations of Ii>i(2430) and Di(2420) 
creation operators) constitutes an important first step to study decays —)• Di + I + v 

using lattice QCD. Such a study is of particular interest, because there is a long standing 
conflict between theory and experiment regarding the corresponding branching ratios (“ 1/2 
versus 3/2 puzzle”). A solid understanding of these decays is in turn necessary for a precise 
determination of the standard model parameter Vcb (cf. [4T] for a detailed discussion). The 
results and corresponding techniques discussed in section 15.1.21 can be used to extend existing 
lattice computations of decays —>■ D** + 1 + 1 /, where D** is currently limited to = O'*" 

and = 2'^, but does not include the two Di states [IH |33] . 


A Computation of correlation matrices using the one-end trick 


To compute the elements of the correlation matrices defined in (1141) we first insert the definition 
of the meson creation operators ([5]), 


= t/(r;r +Ar;t)rj(Ar)(,Sy(2))(r + Ar,t) 

' ^ r Ar=±ea;,±e^,±e2 

j;(xW5)(s, 0) C/(s; s + As; 0)Tk{As){Sx^^'>){s + As, 0)|O), 

s As=±ea; ,±ey,±ez 

where U denote APE smeared gauge links. Gaussian smearing is a linear operation on the quark 
fields and can, therefore, be written in terms of a matrix S in color and position space, where 
S = . After writing the vacuum expectation value (0| ... |0) as a path integral and integrating 

over the quark fields, one obtains an average over gauge link configurations (denoted by (...)), 
which includes quark propagators (/ € {u, d, s"'', s“, c+, c“} is the quark flavor). 


^G;rfc;x(i)x(2)W “ 


1 

V/a^ 


E E E E 

r Ar=±ea;,±ej,,±e2 s As=±ea;,±ej,,itez 


^Ti’spin,color^rfc(As)75C/(s;s + As; 0)(^(Zl^^)) ^’'fS’)(s + As, 0; r + Ar, t)t/(r + Ar; r; t) 
757orj(Ar)7o(S’(T)(^))"^S’)(r,t;s,0)^ (25) 


We have also used twisted mass 75 hermiticity, = 75 (Zl*'-^^)“^’^ 75 , where the bar on top 

of the flavor index indicates a flip of the sign in front of the twisted mass term (e.g. if 2 denotes 
u, 2 denotes d, if 2 denotes s+, 2 denotes s~, etc.). Finally, it is convenient to rearrange the 
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expression and to write spin indices A, B,C,D explicitly, 




~vjcfi ^ (757or](Ar)7o)AS (rfc(As)75)cD 

Ar=±e3;,±ey,±e2 As=±e2;,±e^,±e2 

ee( Trcoiorft^(s;s + As; ^’^S)da{s + As,0;r + Ar,t)U{r + Ar;r;t) 




(5(D«)-i5W(r,t;s,0) 


(26) 


To estimate (j26p stochastically, we generate for each gauge link configuration 4 spin diluted 
stochastic timeslice sources (index B = 1,2,3,4) on a randomly chosen timeslice (here wd.o.g. 
at time t = 0 ), 

= 6t,oSABN<^ir) (27) 

(a is a color index, ^4 is a spin index), where the entries of the noise vector A^“(r) are randomly 
and uniformly chosen numbers ±l/\/2 ± i(y/2. Then we solve the 4 linear systems (index 
C = l,2,3,4) 

y 

-A ()>(^’'^)(r,t) = '^{{D^'^'^)-^S){r,t;x)i^^\x) (28) 


and the 4x6 linear systems (indices D = 1, 2,3,4 and As = ± 63 ,, ±ej^, ±e^) 

'^D^‘^\x-,y)(j)^^’^\y) = ^5(x;y + (0,As))[/(y + (0,As);y)^(^)(y) 
y y 

(/>(^’-^’A®)(r,t) = ^((T>(2))"^5)(r,t;y + (0,As))[/(y + (0,As);y)^(^)(y) 
y 

-A (<)-(2.^.^«))t(r, t) = j;(^(^))t(y)t/(y; y + (0, As))( 5 (Il( 2 ))-i,t)(y + (q, As); r, t) (29) 

y 

with respect to </> using [7S] (to minimize the required computation time, we always choose the 
quark flavors such, that the lighter quark corresponds to flavor 1 and the heavier to quark to 
flavor 2). Using these results we compute the quantities 

/((</.(2’''’^^))^5)^(r + Ar,t)U(r + Ar; r; t)(5</.(i’^))B(r, t)\ (30) 

^ ^ r \ / MC 

((• • •)mc denotes the average over the finite number of Monte Carlo generated gauge link con- 
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figurations [cf. Tabled]), for which one can show 


X 


ABCD;ArAs 
X(l)x(2) 


(*) = 


= EE( Trcoior ( 17(s; s + As; 0) {S^’^S)DAis + As, 0; r + Ar, t) 

r s ' ^ 

U{r + Ar; r; t){S{D^^^)~^S)Bcir, t] s, 0 )^ ^ + noise. 


( 31 ) 


where 

^X“(r)(X'’)*(s)^^^ = + noise (32) 

has been used and “noise” denotes unbiased stochastic noise, which decreases proportional to 
1/y/N, where N is the number of gauge link confignrations used to compute (. ..)mc- This 
technique of stochastic estimation is referred to as one-end trick (cf. e.g. (76] [77]). It is most 
efficient for large spatial volnmes and light quark masses. For the computations of D mesons, 
Dg mesons and charmonium states done in this work it has been found to be snperior compared 
to the traditional technique of using point sources and point-to-all propagators [39] . 

To obtain arbitrary elements of the correlation matrices (jl4l) . one simply has to combine (j26p 
and (l3T]l . 

= (757or](Ar)7o)AS X] (rfc(As) 75 )cL>X^(f)^g’^'’^®(t). 

' Ar=±e3;,±ey,±e2 As=±e2;,±e^,±e2 

(33) 
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